Supersolid Order of Frustrated Hard-Core Bosons in a Triangular Lattice System 
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We numerically demonstrate that a supersolid phase exists in a frustrated hard-core boson system 
on a triangular lattice over a wide range of interaction strength. In the infinite repulsion (Ising) 
limit, we establish a mapping to the same problem with unfrustrated hopping, which connects 
the supersolid to the known results in that case. The weak superfluidity can be destroyed or 
strongly enhanced by a next nearest neighbor hopping term, which provides valuable information 
for experimental realization of a supersolid phase on optical lattice. 
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Introduction: A supersolid phase is a state of mat- 
ter exhibiting both diagonal and off-diagonal long-range 
order (ODLRO) Q. Recent possible observation of a su- 
persolid phase y in ^He under pressure has attracted a 
lot of interest. While the microscopic conditions under 
which clean ^He could be in a supersolid phase are still 
unclear, supersolidity is established for hard-core bosons 
on a triangular lattice, which is the focus of many re- 
cent studies |3i, |4, iSi, 1^, H- The supersolid phase is an 
example of ordering by disorder, demonstrated for hard- 
core boson system with unfrustrated nearest neighbor 
(NN) hopping and strong repulsion 0, [^, Q based on ex- 
tensive Quantum Monte Carlo (QMC) simulations and 
theoretical analysis. These theoretical works are mo- 
tivated in part by experimental realizations of lattice 
bosons in ultra-cold atom traps Intriguingly, the su- 
perfluid density in the supersolid phase is found to be 
very small, possibly indicating that the system is near 
a phase boundary[9| to an insulating phase. It is thus 
highly desirable to examine the stability of the super- 
solid phase in more extended models to determine the 
relevant perturbation and possibly to suggest a route of 
getting into a deep supersolid phase for experiment. 

The nature of the state for the hard-core bosons with 
frustrated NN hopping on triangular lattice is another 
open issue, where the model can be mapped to the spin- 
1/2 XXZ antiferromagetic (AF) Heisenberg model which 
suffers from the sign problem. Historically this model 
was the first candidate proposed to realize a spin liquid 
ground state although it turns out to still exhibit 
a three sub-lattice AF long-range-order (LRO) in gen- 
eral, which may persist to large limit [llj. However, 
extensive numerical studies have been limited to near the 
SU(2) point and the precise nature of the ordering 
at larger Jz (or the strong NN repulsion limit for the cor- 
responding boson model) has not been well understood. 

In this Letter, we present a systematic density-matrix 
renormalization group (DMRG) and exact diagonaliza- 
tion (ED) numerical studies of the half-filled ground state 



of the frustrated model over a wide range of the NN re- 
pulsion. We show that a robust supersolid phase does 
exist in this model, and it can be related to the well- 
known supersolid phase of an unfrustrated NN hopping 
model in infinite repulsion (Ising) limit, where a precise 
mapping between the two models by a sign transforma- 
tion can be established. Furthermore, we reveal that the 
supersolid phase is close to a transition to an insulat- 
ing phase, and correspondingly the superfluidity can be 
strongly enhanced (or easily switched off) by tuning a 
next nearest neighbor (NNN) hopping term. 

Supersolid Phase from Isotropic SU(2) Point to Large 
Jz Limit: We consider a simple model of hard-core 
bosons at half-filling on the triangular lattice interact- 
ing via a repulsive term. 
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where 6t is a boson creation operator and rij a boson 
number operator, (ij) denotes NN sites and we shall 
mainly consider the frustrated boson hopping at t < 0. 
This boson Hamiltonian is equivalent to, by a standard 
mapping from hard-core bosons to S=l/2 spins, an XXZ 
Heisenberg model on the triangular lattice: 
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(SjSj- + h.c.) + JzS^S' 



with Jz = V and Jj_ = —2t. In the spin language, a 
charge density wave (solid) order implies infinite-range 
correlations of the z-component of spins at a nonzero 
wavevector, while the superfluid order is equivalent to the 
in-plane ordering of the spins also at a nonzero wavevec- 
tor for the frustrated system. In the following analysis, 
we will work interchangeably in terms of bosons and spin 
variables. 

We first present the numerical results based on DMRG 
method [Hi for a triangular lattice with the total number 
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FIG. 1: (color online) The order parameter squared, {ruz)^ 
and (mx)^, as functions of the anisotropy A = Jz/J±, are 
shown in (a) and (b) respectively, with A'' = 6x6, 9x6 and 
12 X 9. The insets are the corresponding extrapolations in the 
thermodynamic limit. Examples of finite-size scaling of the 
order parameters are also shown in (cl) and (c2) with system 
size up to A*' = 9 X 18. 



of sites N = Ni x N2- We keep up to to = 4096 states in 
each DMRG block for most systems, and the truncation 
error is of the order or less than 10^^. We make use of 
the periodic boundary condition to reduce the finite-size 
effect for a more reliable extrapolation to the thermo- 
dynamic limit. To analyze the magnetic properties of 
the system, we calculate the structure factors Sz{q) and 
5i.(q) defined by ^,(q) = ^E^.^e-''^^'"'-"'^ (Sf 5/) and 

The obtained Sz{ci) and S'j^(q) show Bragg peaks at 
the corners of the hexagonal Brillouin zone [e.g., at go = 
(±47r/3,0)]. In particular, at small A = Jz/J± < 1, 
the peak of S'z(qo) is very weak, while that of S'_L(qo) 
is very sharp, representing the dominant AF correlation 
in the XY plane. With the increase of A, S'z(qo) grows 
continuously and its value becomes bigger than the in- 
plane ones passing the point A = 1. One can obtain the 
magnetic order parameters based on the finite-size scaling 
of the peak values of S'z(qo) and -S'_L(qo). Specifically, the 
average magnetization (to^) and {mj_) can be determined 
by (to,)^ = SM/N and {m^f = S^{qo)/N, which 
are shown vs. A at = 36, 54 and 108 in the main 
panel of Fig. [TJa) and (b). 

Nonzero (to^)^ and (to_l)^ in the thermodynamic limit 
will correspond to the diagonal LRO and ODLRO, re- 
spectively. Examples of the finite-size scaling are shown 
in FiglH (cl) and (c2) by plotting (m^)^ and {m±)'^ as 



FIG. 2: (color online) Superfluid stiffness ps (in units of J±) 
as a function of the anisotropy A — Jz/J± obtained from ED 
(iV = 36) or DMRG {N = 54) calculations. 



functions of 1/A^. Thus obtained order parameters ex- 
trapolated to the thermodynamic limit are presented in 
the insets of Fig. [T](a) and (b). At small A, the system is 
in the pure superfluid phase with magnetic order solely 
lying in the XF-plane with (toj_)^ 7^ and (to^) = 
(c/. the inset of Fig. [T]). By contrast, at A > 1 the 
thrce-sublattice AFLRO develops in both z-dircction and 
XY plane, characterized by nonzero values of {rriz)'^ and 
(tox)^. Here (to^) monotonically increases with A from 
the isotropic point (A = 1), consistent with the spin- wave 
picture of co-planar ordering in the XZ plane. For the 
corresponding boson system, our results suggest a supcr- 
solid phase with coexisting diagonal LRO and ODLRO 
at A > Ac. The phase boundary Ac between the su- 
perfluid phase and supersolid phase is very close to the 
isotropic point Ac = {Jz/Ji_)c = iV/2\t\)c 1.00. 

We note that at large Jr case, e.g., A = 8, (mz)'^ ~ 
0.049 {{mz) ~ 0.24) which is much larger than the corre- 
sponding value at the isotropic point, while the XK-plane 
magnetization reduces to {m±)'^ = 0.0036 {{m±) = 0.06) 
in the thermodynamic limit. Though this superfluid or- 
dering is small, its value actually is comparable with that 
of the unfrustratcd hard-boson supersolid^, 0, @| at the 
same large Jz limit as we have checked numerically. In 
the following wc can further establish the presence of 
the superfluidity through the calculation of the superfluid 
density ps by adding a nonzero twist phase at the system 
boundary with Ps = ^ oc [E{e^ = tt) - E{e^ = 0)]. 
We obtain ps by adding a twist phase = tt in both 
ED and DMRG calculations which are shown in Fig. 
[5] as a function of A for A^ = 36 and A^ = 54. From the 
figure, we can see that the finite-size effect for ps is very 
weak and ps should remain finite in the thermodynamic 
limit, consistent with the finite {m±)'^ in Fig. [1] At 
Jz > 10, the DMRG becomes difficult to converge as the 
energy from Jz term becomes dominant, the ED results 
further suggest that the obtained ps should be nonzero 
over the whole range of Jz with a value matching with the 
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FIG. 3: (color online) Illustration of the paths contributing 
to the high temperature expansion, and the associated sign 
rule. In (a), we show a representative loop appearing in the 
expansion of the partition function. The ± signs show one of 
two alternating spin configurations allowed around the loop, 
which allows only even-length loops. In (b) , we show a repre- 
sentative path in the expansion of the off-diagonal correlation 
function between sites i and j. For such a path with an odd 
number of steps, the auxiliary sites i + a, j + a (shown in 
circles) must be antiparallel. 



ones for the unfrustrated bosons in the large limit [y]. 

Supersolid Order in the Ising Limit: Now we turn to 
the interesting limit of A — > c», where the direct in-plane 
magnetic ordering in the numerical results is very weak 
and the SW theory suggests that it vanishes as square 
root of J±/Jz- Clearly, here a rigorous examination is 
desired. At — > oo, the XXZ model in Eq. (2) reduces 
to the classical Ising AF on a triangular lattice. This 
classical model is well-known to have a macroscopic de- 
generacy of ground states, which correspond to all spin 
configurations with exactly one frustrated bond per tri- 
angle In this limit, the XXZ model reduces to 
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where Pc is a projection operator onto the classical Ising 
ground state manifold (IGSM). 

The unfrustrated model with Jj^ < has been stud- 
ied previously, and shown to exhibit supersolid order. 
We now show how supcrsolidity in the frustrated case 
Jj^ > can be deduced from those known results. We 
first consider thermodynamic and other properties that 
can be deduced from the partition function and diagonal 
expectation values of the form 
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where O is any function of the z-components of the spins, 
or boson occupation numbers. From such quantities we 
can calculate the free energy, and the diagonal (solid) cor- 
relations. We show that Z is an even function of J±, and 
so these properties are identical for the frustrated and 
unfrustrated cases. To see this, consider the high tem- 
perature expansion of Z in powers of f3J± . The terms in 



the expansion consist of successive actions of bond oper- 
ators of the form PcS^ SJ Pc on nearest-neighbor links, 
with a factor of I3J± accompanying each bond operator. 
To achieve a non-zero expectation value in the trace, the 
boson number on each site must be unchanged after the 
action of all these operators. Graphically, we may rep- 
resent each factor of Sf S~ on the lattice as an arrow 
pointing from site j to site i, and we require this "vector 
field" have zero divergence, i.e. the arrows close into "ex- 
change" loops. Now consider the contribution from any 
particular state in the trace. Due to the projection, each 
bond operator has a non-zero action only if i and j are 
"flippable" , i.e. the two other spins on each triangle con- 
taining i or j are anti-parallel. Now let us circumscribe 
each exchange loop on our graphical representation by a 
neighboring loop as in Fig. [3] In order that all sites on 
the exchange loops are flippable, spins on the neighbor- 
ing loops must alternate, which requires that all of the 
neighboring loops must have an even number of sites. 
This in turn requires that the total number of links on 
each exchange loop is even. Thus Z\0: J\_] is indeed an 
even function of J^_. 

Now consider the off-diagonal correlation function. 
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Once again, one may consider the high temperature ex- 
pansion of the numerator (the denominator Z has al- 
ready been shown to be even). In this case, contribu- 
tions must be divergenceless except at the sites i and j, 
which appear as source and sink, respectively. One can 
understand the behavior by considering just the simplest 
terms, in which the arrows form a single path connecting 
i to j (see Fig.|3]). Now form a tightly circumscribing loop 
about this path. As above, for any state in the trace to 
contribute, the spins on the sites of this neighboring 
loop must alternate. Moreover, the 6 spins neighboring 
i and j must also alternate since Sj'S~ acts directly on 
these states. Now consider the product 4S'f^_^S'|_|_g, where 
d is any nearest-neighbor vector, acting on a state which 
contributes to the trace. Because of the alternating spins 
around the circumscribing loop, this factor gives the par- 
ity of the exchange path, i.e. it equals -1-1 for an even 
path and —1 for an odd path. One may show that this 
conclusion is unaffected by additional closed loops, which 
appear as higher order terms in the high temperature se- 
ries. Since this is true for every term in the expansion, 
we find {S+S-)\j^yo = {ASt^^S^^^^S+ S;)\^.j^. 

The above observations lead us to the conclusion that 
the supersolid phase survives even for the frustrated 
hard-core boson system at strong repulsion (V = Jz ^ 
oo) limit as it maps to the unfrustrated model[l, 0, 0]. 
ED calculation of the projected Hamiltonian gives the 
energy per site varying between — 0.1721< (TV = 36) to 
— 0.1678t {N = 72) with a possible extrapolating value 
E ~ —(0.162 ± 0.005)t in the thermodynamic limit. 
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FIG. 4: (color online) Superfluid stiffness ps (in units of J±) 
vs. J' /J± for systems with NNN hopping t' (superechange 
J'). The finite-size order parameters {niz)'^ and (mx )2 at 
TV = 54 are shown in the inset. 



Enhancement of the Superfluidity and the Ordering of 
the Supersolid Phases: To understand the underlying 
reason why the superfluid stiffness is relatively weakQ, 
we add a NNN hopping t' = —J'/2 term. For simplicity 
we only present the results for = — 2i > and large 
Jz (Ising) limit, although the obtained results also apply 
to both models with a finite range of Jz ■ 

The superfluid stiffness is calculated using the ED 
method for the projected Hamiltonian. As illustrated 
in Fig. m Ps is relatively small at J' = compared to its 
value in the region J' /J± < (negative sign represents 
an unfrustrated NNN hopping). In fact, ps monotoni- 
cally increases with a negative J' and when J' — — J_l, 
Ps becomes comparable to the value of a pure superfluid 
phase (i.e., the ferromagnetic XY model on triangular 
lattice[6|). Clearly a boson system at J' = is indeed 
near the phase boundary of an insulating phase, which 
occurs at {J'/J±)f. ^ 0.2 (which we further identify as 
a solid phase with diagonal LRO). The finite-size order 
parameters (m^)^ and {m±y for = 54 are also shown 
in the inset of the Fig. [H where the enhancement of 
the peaks of the structure factors Sz{qo) and S±{qo) are 
clearly seen as we continuously turn on the negative J'. 
Thus the resulting phase is a supersolid phase with strong 
diagonal LRO and superfiuidity. 

These observations, and the precise nature of the su- 
persolid ordering, can be rationalized by simple energetic 
arguments in the large Jz limit. For the NN hopping 
case (J' = 0), the constraint that neighboring spins to 
the hopping path must alternate tends to enhance hop- 
ping which takes "60°" turns (forming hexagonal path), 
which keeps the bosons on two of the three sublattices. 
Moreover, the third sublattice, on which hopping does 
not proceed, must be substantially polarized. Thus the 
three sublattice ordering, (Sz) =(— m, — m, 2m + S) is 
favored energetically, consistent with a ferrimagnetic or- 
dered phase 0, When the NNN hopping term is 



dominant (—J' > J±), bosons tend to hop on a single 
(say A) sublattice, while spins on the neighboring sites 
from sublattices B and C are individually preferred to be 
ferromagnetically aligned, with B and C spins antiparal- 
lel to each other. This corresponds to (Sz) =(0, m, —to) 
or "antiferromagnetic" ordering, which we therefore ex- 
pect in the large — J' limit. This is indeed supported by 
numerics, which will be presented elsewhere. 

In summary, we have established a robust supersolid 
phase for the frustrated hard-core bosons on triangular 
lattice at half-filling based on extensive numerical cal- 
culations and analytical analysis. The observed super- 
solidity is an example of ordering by disorder elegantly 
realized for such a frustrated system. Furthermore, we 
have found that the supersolid phases for the hard-core 
boson models with only NN hoppings are quite close to 
a pure solid phase in both frustrated and unfrustrated 
cases. But a small unfrustrated NNN hopping term can 
push the boson systems into a deep supersolid phase with 
greatly enhanced superfluidity. Our theoretical study can 
thus provide a solid foundation for the experimental re- 
alization of supersolid phase on triangular optical lattice. 
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